package org.platform.hrv.lib.hrv.calc.parameter.v2;

import lombok.var;
import org.apache.commons.math3.analysis.interpolation.LinearInterpolator;
import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction;
import org.apache.commons.math3.complex.Complex;
import org.apache.commons.math3.stat.StatUtils;
import org.apache.commons.math3.stat.descriptive.DescriptiveStatistics;
import org.apache.commons.math3.transform.DftNormalization;
import org.apache.commons.math3.transform.FastFourierTransformer;
import org.apache.commons.math3.transform.TransformType;
import org.platform.hrv.lib.common.HRVMetrics;
import org.platform.hrv.lib.hrv.calc.parameter.HRVDataProcessor;
import org.platform.hrv.lib.hrv.calc.parameter.HRVParameterEnum;

import java.util.*;
import java.util.stream.Collectors;

/**
 * HRVDataAllProcessor
 *
 * @program: ecg-zoo
 * @description: All HRV parameters are calculated
 * @author: Mr.Nanke
 **/
public class HRVDataAllProcessor {

    /**
     * Calculates the all hrv parameters of the given data.
     *
     * @param data
     *            Data to calculate the all hrv parameters from
     * @return all hrv parameters
     */
    public static HRVMetrics process(List<Integer> rrIntervals) {
//        List<HRVDataProcessor> processors = new ArrayList<>();
//        var HRVMetrics = new HRVMetrics();
        // 封装成json输出
        double[] nn20Results = HRVDataAllProcessor.calculateNN20(rrIntervals);
        double nn20 = nn20Results[0];
        double pnn20 = nn20Results[1];

        double meannni = calculateMeanNNI(rrIntervals);
        double rmssd = calculateRMSSD(rrIntervals);
        double sdnn = calculateSDNN(rrIntervals);

        Map<String, Double> spectral = HRVDataAllProcessor.calculateFrequencyDomainMetrics(rrIntervals);

        return HRVMetrics.builder()
                .uid("")
                .meanNni(meannni)
                .rmssd(rmssd)
                .sdnn(sdnn)
                .cvsd((rmssd / meannni))
                .cvnni((sdnn / meannni) * 100)
                .rangeNni(calculateRangeNNI(rrIntervals))
                .nni50(calculateNN50(rrIntervals))
                .pnni50(calculatePNN50(rrIntervals))
                .nni20(nn20)
                .pnni20(pnn20)
                .tinn(calculateTINN(rrIntervals))
                .hrvLn(Math.log(rmssd))
                .hrvDiff(sdnn)
                .triangularIndex(calculateTriangularIndex(rrIntervals))
                .sdsd(calculateSDSD(rrIntervals))

                .meanHr((60000.0 / meannni) > 200 ? 200 : (60000.0 / meannni))
                .maxHr(calculateMaxHR(rrIntervals))
                .minHr(calculateMinHR(rrIntervals))
                .stdHr(calculateStdHR(rrIntervals))

                .lfnu(spectral.get("lfnu"))
                .hfnu(spectral.get("hfnu"))
                .lfHfRatio(spectral.get("lf_hf_ratio"))
                .lf(spectral.get("LF"))
                .hf(spectral.get("HF"))
                .vlf(spectral.get("VLF"))
                .totalPower(spectral.get("total_power"))

                .build();
    }


    // 1/128秒的毫秒数
    private static final double BIN_WIDTH = 7.8125;

    // 基础时域指标
    public static double calculateMeanNNI(List<Integer> rrIntervals) {
        return rrIntervals.stream().mapToDouble(Integer::intValue).average().orElse(0.0);
    }

    public static double calculateSDNN(List<Integer> rrIntervals) {
        DescriptiveStatistics stats = new DescriptiveStatistics();
        rrIntervals.forEach(stats::addValue);
        return stats.getStandardDeviation();
    }

    public static double calculateRMSSD(List<Integer> rrIntervals) {
        if (rrIntervals.size() < 2) return 0.0;
        double sum = 0;
        for (int i = 1; i < rrIntervals.size(); i++) {
            double diff = rrIntervals.get(i) - rrIntervals.get(i-1);
            sum += diff * diff;
        }
        return Math.sqrt(sum / (rrIntervals.size()-1));
    }

    // 高级时域指标
    public static double calculateCVNNI(List<Integer> rrIntervals) {
        double mean = calculateMeanNNI(rrIntervals);
        double sdnn = calculateSDNN(rrIntervals);
        return (sdnn / mean) * 100;
    }

    public static double calculateRangeNNI(List<Integer> rrIntervals) {
        IntSummaryStatistics stats = rrIntervals.stream()
                .mapToInt(Integer::intValue)
                .summaryStatistics();
        return stats.getMax() - stats.getMin();
    }

    public static int calculateNN50(List<Integer> rrIntervals) {
        int count = 0;
        for (int i = 1; i < rrIntervals.size(); i++) {
            if (Math.abs(rrIntervals.get(i) - rrIntervals.get(i-1)) > 50) {
                count++;
            }
        }
        return count;
    }

    public static double calculatePNN50(List<Integer> rrIntervals) {
        int total = rrIntervals.size() - 1;
        if (total == 0) return 0.0;
        return (calculateNN50(rrIntervals) / (double) total) * 100.0;
    }

    // 直方图分析指标
    public static Map<Double, Integer> buildHistogram(List<Integer> rrIntervals) {
        double min = Collections.min(rrIntervals);
        double max = Collections.max(rrIntervals);
        int numBins = (int) Math.ceil((max - min) / BIN_WIDTH);

        Map<Double, Integer> histogram = new TreeMap<>();
        for (double bin = min; bin <= max; bin += BIN_WIDTH) {
            final double currentBin = bin;
            long count = rrIntervals.stream()
                    .filter(v -> v >= currentBin && v < currentBin + BIN_WIDTH)
                    .count();
            histogram.put(bin, (int) count);
        }
        return histogram;
    }

    public static double calculateTriangularIndex(List<Integer> rrIntervals) {
        Map<Double, Integer> histogram = buildHistogram(rrIntervals);
        int maxCount = Collections.max(histogram.values());
        return (double) rrIntervals.size() / maxCount;
    }

    public static double calculateTINN(List<Integer> rrIntervals) {
        Map<Double, Integer> histogram = buildHistogram(rrIntervals);
        // 寻找最大频数对应的bin
        Map.Entry<Double, Integer> maxEntry = Collections.max(histogram.entrySet(),
                Map.Entry.comparingByValue());

        // 向左找第一个最小值
        Double left = maxEntry.getKey();
        while (histogram.containsKey(left - BIN_WIDTH)) {
            if (histogram.get(left - BIN_WIDTH) > histogram.get(left)) break;
            left -= BIN_WIDTH;
        }

        // 向右找第一个最小值
        Double right = maxEntry.getKey();
        while (histogram.containsKey(right + BIN_WIDTH)) {
            if (histogram.get(right + BIN_WIDTH) > histogram.get(right)) break;
            right += BIN_WIDTH;
        }

        return right - left + BIN_WIDTH;
    }

    // 心率相关指标
    /**
     * 心率均值
     * @param rrIntervals
     * @return
     */
    public static double calculateMeanHR(List<Integer> rrIntervals) {
        return 60000.0 / calculateMeanNNI(rrIntervals);
    }

    public static double calculateMinHR(List<Integer> rrIntervals) {
        double maxRR = Collections.max(rrIntervals);
        return 60000.0 / maxRR;
    }

    public static double calculateMaxHR(List<Integer> rrIntervals) {
        double minRR = Collections.min(rrIntervals);
        return 60000.0 / minRR;
    }

    /**
     * 心率标准差
     * @param rrIntervals
     * @return
     */
    public static double calculateStdHR(List<Integer> rrIntervals) {
        List<Double> heartRates = rrIntervals.stream()
                .map(rr -> 60000.0 / rr)
                .collect(Collectors.toList());
        DescriptiveStatistics stats = new DescriptiveStatistics();
        heartRates.forEach(stats::addValue);
        return stats.getStandardDeviation();
    }

    // 新增FFT长度计算相关方法
    private static int nextPowerOfTwo(int n) {
        if (n <= 0) return 1;  // 关键修复：处理0和负值
        return (int) Math.pow(2, Math.ceil(Math.log(n) / Math.log(2)));
    }


    private static final double SAMPLING_RATE = 4.0; // Hz
    private static final double[] FREQ_BANDS = {0.0033, 0.04, 0.15, 0.4}; // VLF, LF, HF边界

    // 对信号应用汉宁窗
    private static double[] applyHanningWindow(double[] signal) {
        // 创建一个新的数组，用于存储应用汉宁窗后的信号
        double[] windowed = new double[signal.length];
        // 遍历信号数组
        for (int i=0; i<signal.length; i++) {
            // 计算汉宁窗的系数
            double factor = 0.5 * (1 - Math.cos(2*Math.PI*i/(signal.length-1)));
            // 将信号乘以汉宁窗系数
            windowed[i] = signal[i] * factor;
        }
        // 返回应用汉宁窗后的信号
        return windowed;
    }

    // 计算NN20和pNN20
    public static double[] calculateNN20(List<Integer> values) {
        int count = 0;
        for (int i = 1; i < values.size(); i++) {
            if (Math.abs(values.get(i) - values.get(i - 1)) > 20) {
                count++;
            }
        }
        double pnn20 = (count * 100.0) / (values.size() - 1);
        return new double[]{count, pnn20};
    }

    // 计算相邻差值的标准差 (SDSD)
    public static double calculateSDSD(List<Integer> values) {
        double[] diffs = new double[values.size() - 1];
        for (int i = 1; i < values.size(); i++) {
            diffs[i - 1] = values.get(i) - values.get(i - 1);
        }
        double meanDiff = Arrays.stream(diffs).average().orElse(0.0);
        double sum = 0;
        for (double d : diffs) {
            sum += Math.pow(d - meanDiff, 2);
        }
        return Math.sqrt(sum / diffs.length);
    }

    // 其他指标
    public static double calculateCVSD(List<Integer> rrIntervals) {
        double meanNNI = calculateMeanNNI(rrIntervals);
        double rmssd = calculateRMSSD(rrIntervals);
        return (rmssd / meanNNI);
    }

    public static double calculateHRVDiff(List<Integer> rrIntervals) {
        return calculateSDNN(rrIntervals);
    }

    public static double calculateHRVLn(List<Integer> rrIntervals) {
        return Math.log(calculateRMSSD(rrIntervals));
    }


    //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------**

    // 计算频域指标
    public static Map<String, Double> calculateFrequencyDomainMetrics(List<Integer> rriData) {
        Map<String, Double> metrics = new LinkedHashMap<>();

        // 对RRI数据进行重采样，使其等间隔
        double samplingRate = 5.0; // Hz
        double[] resampledData = resampleRRIData(rriData, samplingRate);

        // 去除直流分量（均值）
        double mean = StatUtils.mean(resampledData);
        for (int i = 0; i < resampledData.length; i++) {
            resampledData[i] -= mean;
        }

        // 应用汉宁窗
        double[] windowedData = applyHanningWindow(resampledData);

        // 确保数据长度是2的幂（FFT要求）
        int originalLength = windowedData.length;
        int paddedLength = nextPowerOfTwo(originalLength);
        double[] paddedData = Arrays.copyOf(windowedData, paddedLength);

        // 执行FFT
        FastFourierTransformer transformer = new FastFourierTransformer(DftNormalization.STANDARD);
        Complex[] fftResult = transformer.transform(paddedData, TransformType.FORWARD);

        // 计算功率谱密度
        int n = paddedLength;
        double[] powerSpectrum = new double[n/2];
        for (int i = 0; i < n/2; i++) {
            double real = fftResult[i].getReal();
            double imag = fftResult[i].getImaginary();
            powerSpectrum[i] = (real * real + imag * imag) / n;
        }

        // 计算频带功率
        double frequencyResolution = samplingRate / n;

        // VLF: 0.003-0.04 Hz
        double vlfPower = calculateBandPower(powerSpectrum, frequencyResolution, 0.003, 0.04);

        // LF: 0.04-0.15 Hz
        double lfPower = calculateBandPower(powerSpectrum, frequencyResolution, 0.04, 0.15);

        // HF: 0.15-0.4 Hz
        double hfPower = calculateBandPower(powerSpectrum, frequencyResolution, 0.15, 0.4);

        // Total Power 0.003-0.4 Hz (VLF + LF + HF)
        double totalPower = vlfPower + lfPower + hfPower;

        // 计算归一化功率
        double lfnu = lfPower / (lfPower + hfPower) * 100;
        double hfnu = hfPower / (lfPower + hfPower) * 100;

        // LF/HF比率
        double lfHfRatio = lfPower / hfPower;

        metrics.put("LF", lfPower);
        metrics.put("HF", hfPower);
        metrics.put("lfnu", lfnu);
        metrics.put("hfnu", hfnu);
        metrics.put("lf_hf_ratio", lfHfRatio);
        metrics.put("VLF", vlfPower);
        metrics.put("total_power", totalPower);

        return metrics;
    }


    // 对RRI数据进行重采样
    private static double[] resampleRRIData(List<Integer> rriData, double samplingRate) {
        // 计算累积时间
        double[] timePoints = new double[rriData.size()];
        timePoints[0] = 0;
        for (int i = 1; i < rriData.size(); i++) {
            timePoints[i] = timePoints[i-1] + rriData.get(i);
        }

        // 确定新时间轴
        double totalTime = timePoints[timePoints.length - 1];
        int newLength = (int) (totalTime * samplingRate / 1000.0);
        double[] newTimePoints = new double[newLength];
        for (int i = 0; i < newLength; i++) {
            newTimePoints[i] = i * 1000.0 / samplingRate;
        }

        // 使用线性插值
        LinearInterpolator interpolator = new LinearInterpolator();
        double[] array = rriData.stream().mapToDouble(Integer::doubleValue).toArray();
        PolynomialSplineFunction spline = interpolator.interpolate(timePoints, array);

        double[] resampledData = new double[newLength];
        for (int i = 0; i < newLength; i++) {
            try {
                resampledData[i] = spline.value(newTimePoints[i]);
            } catch (Exception e) {
                // 处理超出范围的情况
                if (newTimePoints[i] < timePoints[0]) {
                    resampledData[i] = rriData.get(0);
                } else {
                    resampledData[i] = rriData.get(rriData.size() - 1);
                }
            }
        }

        return resampledData;
    }

    // 计算频带功率
    private static double calculateBandPower(double[] powerSpectrum, double frequencyResolution,
                                             double lowFreq, double highFreq) {
        int startIndex = (int) (lowFreq / frequencyResolution);
        int endIndex = (int) (highFreq / frequencyResolution);

        // 确保索引在有效范围内
        startIndex = Math.max(0, Math.min(startIndex, powerSpectrum.length - 1));
        endIndex = Math.max(0, Math.min(endIndex, powerSpectrum.length - 1));

        double bandPower = 0;
        for (int i = startIndex; i <= endIndex; i++) {
            bandPower += powerSpectrum[i];
        }

        // 乘以频率分辨率以获得正确的功率值
        bandPower *= frequencyResolution;

        return bandPower;
    }

}